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Abstract 

In this paper we present closed-form solutions for efficiently updating the prin- 
cipal components of a set of n points, when m points are added or deleted from 
the point set. For both operations performed on a discrete point set in W^, we can 
compute the new principal components in 0{m) time for fixed d. This is a signifi- 
cant improvement over the commonly used approach of recomputing the principal 
components from scratch, which takes 0{n + m) time. An important application of 
the above result is the dynamical computation of bounding boxes based on princi- 
pal component analysis. PCA bounding boxes are very often used in many fields, 
among others in computer graphics for collision detection and fast rendering. We 
have implemented and evaluated few algorithms for computing dynamically PCA 
bounding boxes in M^. In addition, we present closed- form solutions for comput- 
ing dynamically principal components of continuous point sets in and M^. In 
both cases, discrete and continuous, to compute the new principal components, no 
additional data structures or storage are needed. 

1 Introduction 

Principal component analysis (PCA) [15] is probably the oldest and best known of the 
techniques of multivariate analysis. The central idea and motivation of PCA is to reduce 
the dimensionality of a point set by identifying the most significant directions (principal 
components). Let P = {pi,P2, ■ ■ ■ ,Pn} be a set of vectors (points) in M'^, and fl = 
{fii, fi2, ■ ■ ■ , fJ'd) £ 1^*^ be the center of gravity of P. For 1 < k < d, we use pi^k to denote 
the k-th coordinate of the vector pi. Given two vectors u and v, we use {u, v) to denote 
their inner product. For any unit vector v G W^, the variance of P in direction v is 
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var(P, v) = - ^{Pt - fl, v)"^. (i; 
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The most significant direction corresponds to the unit vector vi such that var(P, vi) 
is maximum. In general, after identifying the j most significant directions Vi,...,Vj, 
the (j + l)-th most significant direction corresponds to the unit vector vj+i such that 
var(P, Vj+i) is maximum among all unit vectors perpendicular to ?Ti, {72, • • • , vj. 

It can be verified that for any unit vector v G M*^, 



where S is the covariance matrix of P. S is a symmetric d x d matrix where the (z, j)-th 
component, (Jij,l < i, j < d, is defined as 



The procedure of finding the most significant directions, in the sense mentioned above, 
can be formulated as an eigenvalue problem. If Ai > A2 > ■ ■ ■ > are the eigenvalues of 
E, then the unit eigenvector vj for \j is the j-th most significant direction. All AjS are 
non-negative and Xj = vai{X,Vj). Since the matrix E is symmetric positive definite, its 
eigenvectors are orthogonal. 

Computation of the eigenvalues, when d is not very large, can be done in 0{d^) time, 
for example with the Jacobi or the QR method [18]. For very large d, the problem of 
computing eigenvalues is non-trivial. In practice, the above mentioned methods for com- 
puting eigenvalues converge rapidly. In theory, it is unclear how to bound the running 
time combinatorially and how to compute the eigenvalues in decreasing order. In [10] a 
modification of the Power method [IT] is presented, which can give a guaranteed approx- 
imation of the eigenvalues with high probability. 

Examples of many applications of PCA include data compression, exploratory data 
analysis, visualization, image processing, pattern and image recognition, time series pre- 
diction, detecting perfect and reflective symmetry, and dimension detection. The thorough 
overview over PCA's applications can be found for example in the textbooks [13] and [T5] . 
Most of the applications of PCA are non-geometric in their nature. However, there are 
also few purely geometric applications that are quite spread in computer graphics. Ex- 
ample are the estimation of the undirected normals of the point sets or computing PCA 
bounding boxes (bounding boxes determined by the principal components of the point 
set). 

Dynamic versions of the above applications, i.e., when the point set (population) 
changes, are of big importance and interest. In this paper we present closed-form solutions 
for efficiently updating the principal components of a dynamic point set. We also consider 
the computation of the dynamic PCA bounding boxes - a very important application in 
many fields including computer graphics, where the PCA boxes are used to maintain 
hierarchical data structures for fast rendering of a scene or for collision detection. 

Based on the theoretical results in this paper, we have implemented several algorithms 
for computing PCA bounding boxes dynamically. 

The organization and the main results of the paper are as follows: In Section [2] we 
present closed-form solutions for efficiently updating the principal components of a set of 
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n points, when m points are added or deleted from the point set. For both operations 
performed on a discrete point set in M'^, we can compute the new principal components 
in 0{m) time for fixed d. This is a significant improvement over the commonly used ap- 
proach of recomputing the principal components from scratch, which takes 0{n+m) time. 
In Section [3] we consider solutions for the static and dynamic versions of the bounding 
box problem. In Section H] we present and verify the correctness of the theoretical results 
presented in the Chapter [2l We have implemented several dynamic PCA bounding box al- 
gorithms and evaluated their performances. Conclusion and open problems are presented 
in Section [51 In the appendix we consider the computation of the principal components 
of a dynamic continuous point set. We give closed form-solutions when the point set is 
a convex polytope or a boundary of a convex polytope in or M?. When the point 
set is a boundary of a convex polytope, we can update the new principal components 
in 0{k) time, for both deletion and addition, under the assumption that we know the k 
facets in which the polytope changes. Under the same assumption, when the point set 
is a convex polytope in or M'^, we can update the principal components in 0{k) time 
after adding points. But, to update the principal components after deleting points from 
a convex polytope in or we need 0{n) time. This is due to the fact that after a 
deletion the center of gravity of the old convex hull (polyhedron) could lie outside the 
new convex hull, and therefore, a retetrahedralization is needed (see Subsection 16.1.11 and 
Subsection 16.2.11 for details) . 

2 Updating the principal components efficiently - 
discrete case in 

In this subsection, we consider the problem of updating the covariance matrix S of a 
discrete point set P = {'Pi,p2, ■ ■ ■ ,Pn} in K*^, when m points are added or deleted from 
P. We give closed-form solutions for computing the components of the new covariance 
matrix S'. Those closed-form solutions are based on the already computed components 
of S. The main result of this section is given in the following theorem. 

Theorem 1 Let P be a set of n points in Mf^ with known covariance matrix S. Let P' 
he a point set in W^, obtained by adding or deleting m points from P. The principal 
components of P' can be computed in 0{m) time for fixed d. 

Proof. 

Adding points 

Let Pm = {pn+i,Pn+2, • • • , Pn+m} be a poiut Set with center of gravity /I"" = (/i™, /i^, . . . , /i™). 
We add Pm to P obtaining new point set P'. The j-th component, fi'j, I < j < d, of the 
center of gravity fl' = {fi'i, /ig, . . . , /i^) of P' is 
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The (i, j)-th component, cr^ -, 1 < i, j < d, of the covariance matrix S' of P' is 
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Plugging- in the values of /i^ and /i^- in (jl]), we obtain: 
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Since Efc=i(PA:,i - /ij) = 0, 1 < z < d, we have 
Plugging- in the values of fi^ and /x^- in (jSj), we obtain: 
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Since Yltn+iiPKi - /i^) = 0, 1 < z < ci, we have 
where 



^ n+m 
k=n+l 

is the j-th element of the covariance matrix of the point set Pm- 
Finally, we have 

1 TlTtl 

4 = 4,1 + 4,2 = ;^(-^^. + -<) + - - /^D- (8) 

Note that crj^, and therefore a[p can be computed in 0{m) time. Thus, for a fixed 
dimension (i, the covariance matrix E also can be computed in 0{m) time. 

Deleting points 

Let Pm = {Pn~m+i,Pn-m, ■ ■ ■ yPn} be a subset of the point set P, and let /I™ = 
{/j,^ , fi^ , . . . , fi^) be the center of gravity of Pm- We subtract Pm. from P, obtain- 
ing new point set P'. The j-th component, fi'j, 1 < j < d, of the center of gravity 
/^' = if^'i, ■ ■ ■ , f^'d) ofP' is 

/ 1 sr^n—m 

~ (X]fc=l ^ SA:=n-m+l^'fc,i) 

n , . m m 

n—mi^l n—mi^j 



n—m^3 

The (i, j)-th component, cr^ -, 1 < i, j < d, of the covariance matrix S' of P' is 
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Plugging-in the values of /i^ and jjl^ in ([H]), we obtain: 
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Since YTk=AVk,i - /ij) = 0, 1 < z < d, we have 
Plugging-in the values of /i^ and /i^- in (fTOj) . we obtain: 
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Since ELn-m+ib^,* - /ir) = 0, 1 < z < ci, we have 

_/ n I n'^m ( ,,m\( ,, , ,m\ 

^ii,2 = ^(^ij + (^^if^i - f^i AH - /^j (12) 

where 
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is the i, j-th element of the covariance matrix of the point set Pm- 
Finally, we have 

1 TlTfi 

4 = 4,1 + 4,2 = i^^io - "^O - 7 (^i - (/^i - /^D • (13) 

Note that crj^, and therefore a[p can be computed in 0{m) time. Thus, for a fixed 
dimension the covariance matrix E also can be computed in 0{m) time. □ 

As a corollary of ([S]), in the case when only one point, Pn+i-, is added to a point set P, 
the elements of the new covariance matrix are given by 

Tl TL 

^'ij = + 4,2 = -—T^ij + r , 1N2 (P"+1-^ ~ f^i){Pn+l,j - ^i), (14) 

n + i (n + ij^ 
6 



and can be computed in 0(1) time. 

Similarly, as a corollary of fll3l) . in the case when only one point, p^-, is deleted from a 
point set P, the elements of the new covariance matrix are given by 

and also can be computed in 0(1) time. 

The principal components of discrete point sets can be strongly influenced by point 
clusters [12]. To avoid the influence of the distribution of the point set, often continuous 
sets, especially the convex hull of a point set is considered, which lead to so-called con- 
tinuous PCA. Computing PCA bounding boxes [H], [TT], or retrieval of 3D-objects [23], 
are typical applications where continuous PCA are of interest. Due to better readability 
and compactness of the paper, we present the closed-form solutions for dynamic version 
of continuous PCA in the appendix. There, we consider cases when the point set is a 
convex polytope or the boundary of a convex polytope in and M^. 

3 An application - computing PCA bounding boxes 

PCA is a well-established technique for dimensionality reduction and multivariate analy- 
sis, with numerous applications in both static and dynamic context. Examples of many 
applications of PCA include data compression, exploratory data analysis, visualization, 
image processing, pattern and image recognition, time series prediction, detecting per- 
fect and reflective symmetry, and dimension detection (see textbooks [13] and [15] for 
thorough overview over PCA's applications). 

In the rest of this section, we consider solutions for the static and dynamic versions of 
the bounding box problem. 

3.1 Computing bounding boxes - static version 

Substituting sets of points or complex geometric shapes with their bounding boxes is 
motivated by many applications. For example, in computer graphics, it is used to main- 
tain hierarchical data structures for fast rendering of a scene or for collision detection. 
Additional applications include those in shape analysis and shape simplification, or in 
statistics, for storing and performing range-search queries on a large database of samples. 

Computing a minimum-area bounding rectangle of a set of n points in can be 
done in 0(?7.1og?7,) time, for example with the rotating calipers algorithm [22]. O'Rourke 
[T9] presented a deterministic algorithm, a rotating calipers variant in M^, for comput- 
ing the minimum- volume bounding box of a set of n points in . His algorithm re- 
quires 0{n^) time and 0(n) space. Barequet and Har-Peled [1] have contributed two 
algorithms with nearly linear complexity, based on a core-set approach, that compute 
(H-£:)-approximations of the minimum-volume bounding box of point sets in R'^. The 
running times of their algorithms are 0(n + and 0(nlogn -|- n/e^), respectively. 

A further improvement to 0{n + l/£^) running time can be obtained by using a coreset 
of size 0(l/e) by Agarwal, Har-Peled, and Varadarajan [2], and Chan [8]. 
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Numerous heuristics have been proposed for computing a box that encloses a given 
set of points. The simplest heuristic is naturally to compute the axis-aligned bounding 
box of the point set. Two-dimensional variants of this heuristic include the well-known 
R-tree, the packed R-tree [20], the R*-tree [5], the R'^-tree [21], etc. 

A frequently used heuristic for computing a bounding box of a set of points is based 
on PCA. The principal components of the point set define the axes of the bounding box. 
Once the directions of the axes are given, the dimension of the bounding box is easily 
found by the extreme values of the projection of the points on the corresponding axis. 




Two distinguished applications of this heuristic are the OBB-tree [H] and the BOX- 
TREE [3], hierarchical bounding box structures, that support efficient collision detection 
and ray tracing. Computing a bounding box of a set of points in and M.^ by PCA is 
simple and requires linear time. To avoid the influence of the distribution of the point 
set on the directions of the PCs, a possible approach is to consider the convex hull, or 
the boundary of the convex hull CH{P) of the point set P. Thus, the complexity of the 
algorithm increases to 0(n log n). The popularity of this heuristic, besides its speed, lies 
in its easy implementation and in the fact that usually PCA bounding boxes are tight- 
fitting, see Figure [1] for an illustration. Experimental results of the quality of the PCA 
bounding boxes can be found in [II],[I6], and theoretical results in [T2] . 

3.2 Computing bounding boxes - dynamic version 

Dynamic (l-l-e)-approximation bounding box can be efficiently solved by the dynamic data 
structure [9] that can maintain an e-coreset of n points in O(logn) time for any constant 
e > and any constant dimension. That is an improvement of the previous method 
by Agarwal, Har-Peled, and Varadarajan [2] that requires polylogarithmic update time. 
However, both results are more of theoretical importance, since there realization involves 
sophisticated data structures difficult for implementation. 
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Computing PCA bounding boxes of a point set consists of two steps: 1) computing 
the principal components, that define the axes of the bounding box, and 2) computing the 
extremal point along the axes, that determine the size of the bounding box. In Section [2l 
we have presented closed-form solution for efficient update of the principal components 
when we add or delete several points from the point set. In sequel, we consider step 2. 

3.2.1 Computing extremal points 

A trivial way to compute the extremal points is to scan all points, which takes linear time. 

Faster algorithms that compute extremal points dynamically are related to computing 
and maintaining the convex hull of the point set dynamically. Then, one can perform 
extreme point queries in polylogarithmic time. Here, we give an overview of the known 
results of dynamical computation of convex hull and some of related operation on it in 

and M^. 

Brodal and Jacob [B] present a data structure that maintains a finite set of n points 
in the plane under point insertions and point deletions in amortized 0(log?7.) time per 
operation. This data structure requires 0{n) space, and supports extreme point queries in 
a given direction, tangent queries through a given point, and queries for the neighboring 
points on the convex hull in O(logn) time. T. Chan [7] presents a fully dynamic random- 
ized data structure that can answer queries about the convex hull of a set of n points in 
three dimensions, where insertions take O(log^n) expected amortized time, deletions take 
O(log^n) expected amortized time, and extreme-point queries take O(log^n) worst-case 
time. This is the first method that guarantees polylogarithmic update and query cost 
for arbitrary sequences of insertions and deletions, and improves the previous 0(n^)-time 
method by Agarwal and Matousek [T]. 

There are two disadvantages in the above approaches. First, they require a computa- 
tion of the convex hull, which increases the complexity of the algorithms to O(nlogn). 
However, the convex hull computation and building corresponding data structures can be 
done in preprocessing, which is not critical for many applications. Second, the above date 
structures for dynamic convex hull computation are of theoretical importance, they are 
quite difficult for implementation, and to best of our knowledge, they have never been 
implemented. Therefore, in the next section, we consider two simple approaches for com- 
puting extremal points, one is the already mentioned linear scan of all points, and the 
other is a grid approach, refined with several variants. 

4 Practical variants of dynamical PCA bounding boxes 
and experimental results 

The main focus in this section is to show the advantages of the theoretical results presented 
in this paper in the context of computing dynamic PCA bounding boxes. We present three 
practical simple algorithms, and compare their performances. A thorough comparison 
study of different variants of statical PCA bounding boxes the interested reader could 
find in [11]. The algorithms were implemented in C#, C++ and OpenGL, and tested 
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on a Core Duo 2. 33 GHz with 2GB memory. The principal components of all algorithms 
are computed with the closed-form solutions from Section [2l They differ only how the 
extremal points along the principal components are found. The implemented algorithms 
are the following: 

• PCA-AP (PCA-all-points) - finds the extremal points by going through all points. 

• PCA-AGP (PCA-all-grid-points) -the space is discretized by a regular three di- 
mensional axis- aligned grid, with a cube of size e x e x e a.s primer component. See 
Figure [2] for an illustration. The grid size is chosen relatively to the size of the ob- 
ject. Each object is scaled such that its diameter is 1. The values of e are between 
0.001 and 1. The corners of non-empty cells are considered to find the extremal 
points along the principal directions. 




Figure 2: (a) A real world object and its corresponding grid for e = 0.03. Only the 
non-empty grids are visualized, (b) The bounding box of the object obtained by the 
PCA-AGP algorithm. 

• PCA-EGP (PCA-extremal-grid-points) - this is an improvement of the PCA-AGP 
algorithm. To each vertical grid line, i.e., orthogonal to XY plane, two extremal 
corners of the non-empty cell are computed. Thus, we reduced the candidates for 
extremal points from O(^) to O(p-). 

We further reduce the number of points considered in the PCA-AGP and PCA-EGP 
algorithms by replacing the cell corners with the centers of gravity of the cells. Afterwords, 
we expand the resulting box by \^e/2 to ensure that the box contains all original points. 
We have implemented also these variants, but, since for a reasonable big grid size {e > 
0.01) the running time improvements are negligible, we report here only the results of the 
base variants of the algorithms PCA-AGP and PCA-EGP. However, for very dense grid 
the improved version of the both algorithms give better results. 
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In the following experiments, we add (delete) random points from the point set, and 
compare the results of a dynamical versions of PCA bounding boxes with the their corre- 
sponding statical versions (when the covariance matrix of the point set is computed from 
scratch). The time of computing, the volume of a bounding box, and the grid density are 
parameters of interest in this evaluation study. The test were performed on big number of 
real graphics models taken from various publicly available sources (Stanford 3D scanning 
repository, 3D Cafe). Typical samples of the results are given in Table [H Table [21 and 
Table El 

Table 1: Time needed by the PCA bounding box algorithms for the lion model 
(183408 points). The values in the table are the average of results of 100 runs of the 
algorithms, each time adding/ deleting the corresponding number of points. 



Adding/deleting points, e = 0.005 




Ipnt 


Ipnt 


100 puts 


100 puts 


1000 puts 


1000 puts 


algorithm 


static 


dynamic 


static 


dynamic 


static 


dynamic 


PCA-AP 


0.166 s 


0.014 s 


0.171 s 


0.015 s 


0.172 s 


0.016 s 


PCA-AGP 


0.092 s 


0.0095 s 


0.093 s 


0.0085 s 


0.99 s 


0.017 s 


PCA-EGP 


0.0805 s 


0.0055 s 


0.082 s 


0.006 s 


0.092 s 


0.0135 s 



The main conclusions of the experiments are as follows: 

• As expected from the theoretical results, the dynamic versions of the algorithms are 
significantly faster than their static counterparts. Typically, the dynamic versions 
are about an oder of magnitude faster (see Tabled]). 

• The dynamic PCA-AP algorithm is not only significantly faster than its statical 
version, it is also faster than the static version of the PCA-AGP and PCA-EGP 
algorithms. This is due to fact that the brute force manner of finding the extremal 
points is faster than computing the covariance matrix of the new point set from 
scratch, although both algorithms require 0{n) time in the asymptotic analysis. 

• Clearly, the PCA-AGP and PCA-EGP algorithms, that exploit the grid subdivision 
structure, are faster than the PCA-AP algorithm. The price that must be paid for 
this is two- folded. First, an extra preprocessing time for building the grid is needed. 
For the example considered in Tabled], computing the grid takes about 0.4 seconds 
for the PCA-AGP algorithm, and about 0.43 for the PCA-EGP algorithm. Second, 
the resulting bounding boxes are less precise (see Table [2]). 

• As it is shown in Table [3l for grids that are not very sparse {e < 0.03), the approxi- 
mated PCA bounding boxes computed by the PCA-AGP and PCA-EGP algorithms 
are quite close to the exact PCA bounding boxes. 

Tight bounding boxes for the PCA-AGP and PCA-EGP algorithms can be obtained by 
the following approach. Let Pi be the supporting plane at the extremal grid point along 
one principal direction, and let P2 be the plane parallel to Pi, such that the distance 



11 



Table 2: Volume of the PCA bounding box algorithms for the lion model. The values in 
the table are the average of results of 100 runs of the algorithms, each time adding the 
corresponding number of points. 



Adding points, dynamic version, e = 0.005 


algorithm 


Ipnt 


lOpnt 


100 puts 


1000 puts 


10000 puts 


PCA-AP 


285.5 


644.6 


856.3 


1149.1 


1236.4 


PCA-AGP, PCA-EGP 


295.5 


662.7 


880.3 


1221.8 


1263.2 



Table 3: Volumes of the PCA bounding boxes algorithms for lion model for different grid 
density. The values in the table are the average of results of 100 runs of the algorithms, 
each time adding the corresponding number of points. 



Adding 100 points, dynamic version 


algorithm 


e = 0.005 


e = 0.01 


e = 0.03 


e = 0.05 


e = 0.1 


e = 0.2 


PCA-AP 


856.3 


856.3 


856.3 


856.3 


856.3 


856.3 


PCA-AGP, PCA-EGP 


880.3 


904.3 


942.3 


1080.1 


1292.7 


2324.8 



between Pi and P2 is v^e/2, and P2 intersect or is tangent to the grid. We denote by 
S the subspace between Pi and P2. Then, the candidates points for the chosen principal 
direction, that determine the tight bounding box, are all original points that belong to 
cells that have intersection with S. See Fig. [3] for an illustration. However, in the worst 
case all original points have to be checked. 
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Figure 3: For the principal direction PCi, the algorithms PCA-AGP and PCA-EGP 
detect the point gi as extremal grid point, and the point x as extremal point of the 
original point set. However, there are other points (the violet colored circles) that are 
further than x along PCi. 
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Further (theoretical) improvement of the algorithms presented here could be obtained 
if, instead of the point set, we consider its convex hull when we look for extremal points. 
This only makes sense if the convex hull is computed dynamically. Otherwise, computing 
the static convex hull of the points will be more expensive than finding the exact extremal 
points by scanning all points. 




Figure 4: Left: two objects with their PGA bounding boxes. Right: the common PGA 
bounding box. Gomputing the common PGA bounding box dynamically takes 0.004 
seconds, while the static version takes 0.02 seconds. 

4.0.2 Computing efRciently a bounding box of several objects 

An interesting application of the closed-form solutions from Section [2] is to compute the 
principal components of two or more objects with already known covariance matrices. 
Since aij and a™ in ([8]) and f|T5l) are previously known, ct-^- can be computed in 0(1) time. 
Thus, for fixed d the new covariance matrix S and the new principal components can be 
computed also in 0(1) time. This is a significant improvement over the commonly used 
approach to compute the principal components from scratch, which take time linear in 
the number of points. Efficient computation of the common PGA bounding box of several 
object is straightforward. See Fig. |4]for an illustration in M.^. 
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5 Conclusion and future work 



The main contribution of this paper are the closed-form solutions for updating the princi- 
pal components of a dynamic point set. The new principal components can be computed 
in constant time, when a constant number of points are added or deleted from the point 
set. This is a significant improvement of the commonly used approach, when the new prin- 
cipal components are computed from scratch, which takes linear time. The advantages of 
the theoretical results were verified and presented in the context of computing dynamic 
PCA bounding boxes, a very important application in many fields including computer 
graphics, where the PCA boxes are used to maintain hierarchical data structures for fast 
rendering of a scene or for collision detection. We have presented three practical simple 
algorithms and compare their performances. 

In the appendix we consider the computation of the principal components of a dy- 
namic continuous point set. We give closed form-solutions when the point set is a convex 
polytope or the boundary of a convex polytope in or R^. 

An interesting open problem is to find a closed-form solution for dynamical point 
sets different from convex polyhedra, for example, implicit surfaces or B-splines. An 
implementation of computing principal components in a dynamic and continuous setting 
is planned for future work. Applications of the results presented here in other fields, like 
computer vision or visualization, are of high interest. 

There are several further improvements and open problems regarding computing dy- 
namic PCA bounding boxes. Instead of subdividing the space by a simple regular grid, 
one can use more sophisticated data structures, like octrees or binary space partition-trees 
to speed up the time needed to find the extremal points along the principal directions. A 
practical, implementable algorithm for computing the dynamic convex hull of the point set 
(computing extremal point dynamically) would also improve the dynamic PCA bounding 
box algorithms. Finding coresets for dynamic PCA bounding boxes will lead to efficient 
approximation algorithms for PCA bounding boxes. We are also not aware of data struc- 
tures for efficient computation of extremal points both approximately and dynamically. 
Such data structures are also of interest. 
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6 Appendix 

Updating the principal components efficiently - continuous 
case 

Here, we consider the computation of the principal components of a dynamic continuous 
point set. We present a closed form- solutions when the point set is a convex polytope or a 
boundary of a convex polytope in or . When the point set is a boundary of a convex 
polytope, we can update the new principal components in 0{k) time, for both deletion 
and addition, under the assumption that we know the k facets in which the polytope 
changes. Under the same assumption, when the point set is a convex polytope in or 
M^, we can update the principal components in 0{k) time after adding points. But, to 
update the principal components after deleting points from a convex polytope in or 
we need 0{n) time. This is due to the fact that, after a deletion the center of gravity 
of the old convex hull (polyhedron) could lie outside the new convex hull, and therefore, 
a retetrahedralization is needed (see Subsection 16.1.11 and Subsection 16.2.11 for details). 

6.1 Continuous PC A in 

6.1.1 Continuous PCA over a (convex) polyhedron in 

Let P be a point set in M^, and let X be its convex hull. We assume that the boundary 
of X is triangulated (if it is not, we can triangulate it in preprocessing). We choose an 
arbitrary point o in the interior of X, for example, we can choose that o is the center of 
gravity of the boundary of X. Each triangle from the boundary together with o forms a 
tetrahedron. Let the number of such formed tetrahedra be n. The /c-th tetrahedron, with 
vertices xi^k,X2,k,X3^k,X4^k = o, can be represented in a parametric form by Qi{s,t,u) = 
Xi^i + s {xi^i — Xi^i) + 1 (x2,i — X4,i) + u {xs^i — X4^i), for < s, t, u < 1 , and s + t + u < 1 . For 
1 < z < 3, we use Xij^k to denote the i-th coordinate of the vertex xj of the polyhedron 

Qk- 

The center of gravity of the k-th tetrahedron is 

-, _ Jp ~° Jo^~°~' p{Qk{s,t))Qi{s,t) dudtds 

f^fo'-' If 'p(Qk{s,t)) dudtds ' 
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where p{Qk{s, t)) is a mass density at a point Qk{s, t). Since, we can assume p{Qk{s, t)) = 
1, we have 

Tj _ Jo fo'" /o^"°~' Qkis,t)dudtds _ Xi^^.+X2,k+S3,k+X4,k 

^'^ ~ dudtds - 4 

The contribution of each tetrahedron to the center of gravity of X is proportional to its 
volume. If Mfc is the 3x3 matrix whose /-th row is xi^k ~ X4^k, for / = 1 . . . 3, then the 
volume of the k-th tetrahedron is 

\det{Mk)\ 

Vk = volume (Qfe) = 



3! 

We introduce a weight to each tetrahedron that is proportional with its volume, define as 

U!k = ^ = — , 

where v is the volume of X. Then, the center of gravity of X is 

n 

fi=^ WkPk- 
k=l 

The covariance matrix of the fc-th tetrahedron is 

„ _ Sq /o^~° /o^^°~* (Qk{s,t,u)-fl) {Qk{s,t,u)-fl)'^ dudtds 

— (-1 rl-s fl- 



/o So ' lo ' <^'^ '^^ 

^(Ei=i ELife - Mxh,k - pf+ 



The (i, j)-th element of S^, i,j G {1,2,3}, is 

^V,k = ^ ( ELi Eh=l(^i,«,fc ~ f^i){XjXk - + 
YA=li^i,l,k — fJ'i){Xj,l,k — fJ'j)), 



with fl = (/ii, /i2, /is). Finally, the covariance matrix of X is 

^ = Er=i ^i^i' 

with (z, j)-th element 

f^ij = ^ ( ELi Eti Et=i vJi{xi,i,k - f^i)(xj,h,k - f^j) + 

ELl ELi Wii^i,l,k - f^r){Xj,l,k - Pj)^ ■ 

We would like to note that the above expressions hold also for any non-convex poly- 
hedron that can be tetrahedralized. A star-shaped object, where o is the kernel of the 
object, is such example. 
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Adding points 

We add points to P, obtaining a new point set P' . Let X' be the convex hull of P' . 
We consider that X' is obtained from X by deleting n^, and adding tetrahedra. Let 

n Ua rid na riii 

V' = '^Vk + '^Vk-'^Vk = V + '^Vk-'^Vk■ 
k=l k=l k=l k=l k=l 

The center of gravity of X' is 

= Efc=i ^fc/^fc + Efc=i ^fc/^fc - Efc=i ^fc/^fc 

= ^ (ELi ^fc/^fc + Efeli ^fc/^fc - Efc=i ^fc/^fc) (16) 

= ^{vfi+ Ylli Vkfik - Efcii Vkfik) ■ 



Let 



/Ua = — Vkjlk, and = — ^^fc/^fc- 

k=l k=l 

Then, we can rewrite ffT6l) as 



/i = — /i+ /ia - /id- (17) 

The 2-th component of and /Irf, 1 < i < 3, is denoted by fii a and fii^d, respectively. 
The (i, j)-th component, a'^p 1 < 2, j < 3, of the covariance matrix S' of X' is 



4 



^ ( ELi Eti ELi Ki^iAk - f^'i)i^jAk - 
ELi Eti Ki^i,i,k - f^'i)i^j,i,k - + 

^ ( Efcli Eti ELi - f^'i}ixj,h,k - fJ^'j) + 

Y2=i Eti K{xi,uk - f^'i)ixj,i,k - /ij)- 

Efcii Eti ELi ^fc(a;i,i,fc - f^'i){xj,h,k -f^'j)- 

Efc=i Eti ^fc(2;i,«,fc - ^'i){xj,i,k - /i'; 



Let 



where, 



1 



n 4: A 
k=l 1=1 h=l 
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n 4 



^ii,12 = ^^Wk{Xi,l,k - fJ''i){Xj,l,k - ^I'j), (19) 
k=l 1=1 

na 4 4 

^ii,21 = ^^Y1 ^'ki^i'l'k - fJ''i)iXj,h,k - fi'j), (20) 

k=l 1=1 h=l 

ria 4 

f^ii,22 = XI XI ^ ^^'i}i^j,l,k - /^i), (21) 
fc=l 1=1 

= X X X ~ l^'i)(^jAk - fJ^'j), (22) 
k=i 1=1 h=i 

rid 4 

'^ij,32 = X X ~ /^D(a;j,i,fc - /^i)- (23) 

fc=i ;=i 

Plugging-in the values of /i^ and fi'j in f|T8|) . we obtain: 

'^ij,U — J2k=l X]i=l S/i=l '^'k{^i,l,k — ^fJ'i — f^i,a + f^i,d){Xj^h,k — ^fJ'j — l^j,a + ^'j,d) 
= ELl Ef=l Et=l ^fc(a;i,Z,fc - /^i + /^i(l - ^) - /^*,a + ^^^4) 

{Xj,h,k - /ii + /ij(l - ^) - /ij.a + /ii,d) 

= ELl Ef=i EJ=i K{xi,i,k - f^i){xj,h,k - /ij)+ 

ELl Etl ELi KiXi,l,k - /^i)(/ii(l - ^) - /^j,a + /ii,d) + 

ELl Eti ELi w'kif^ii^-v)- + t^i,d)i^jAk - /^i)+ 

ELl ELi ELi Kif^ii^-V)- + /^i,d)(/^i(l - ^) - f^j,a + f^j,d)- 

(24) 

Since X^Li ELi ^'ki^i,i,k - /^i) = 0, 1 < i < 3, we have 

f^'i,ii = V ELl Eti ELi Mxi,i,k - ^^i){xj,h,k - /ij)+ 

^ ELl Etl ELi ^fc(/^i(l - ^) - /^i,a + /^i,d)(/^j(l - ^) - /^j,a + /^i,d) 

= V ELl Etl ELi ^fc(a;i,i,fc - f^i){xj,h,k - /ii) + 

16^(/ii(l - ^) - /ii,a + /ii.d)(/ij(l - :7) - /ii,a + /ij,d)- 

(25) 

Plugging-in the values of fj,[ and fi'j in (fT9|) . we obtain: 
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(^'ij,l2 = ELl E?=l K iHl,k - ^/^i - f^i,a + f^i,d) {Xj,h,k - ^/ij - /ii,a + f^j,d) 

= ELiEti<K;,fc -/^^(a^iAfc (26) 
Since XlLi ELi ^'ki^i,i,k -/ii) = 0,l<i<3, we have 

'^ii,12 = ^ Efc=l ELi '"k{Xi,l,k - fJ'i){Xj,h,k - f^j) + 

^ ELl Eti Mf^ii'^ - ^) - /^*.« + - ^) - + 

= V ELl Eti ^fc(a;i,«,fc - f^i){xj,h,k - 

From (!26l) and (!27|) . we obtain 



(27) 



(28) 

= aij + 20^(/ii(l - ^) - yUi,a + /ii,d)(Aij(l - :^) - /^i.a + /ii,d)- 

Note that cr^ ^ can be computed in 0(1) time. The components cr'^j 21 and cr'-j 22 '^^^ 
be computed in 0{na) time, while 0{nd) time is needed for computing cr^.gi and cr[j .^2- 
Thus, /i' and 

^'ij ~ 'h^^'ij,^^ ~'~ ^ii,12 + '^ii,21 + ^«'i,22 + ^ij,31 + '^ii,32) 

= M(^ii + ^*W+^ii,22 + ^ii,31 + ^^i.32)+ (29) 
^(/^i(l - ^) - /^»,a + /^»,d)(/^i(l - ^) - /^J> + /^J,d) 

can be computed in 0{na + Ud) time. 
Deleting points 

Let the new convex hull be obtained by deleting tetrahedra from and added Ua 
tetrahedra to the old convex hull. If the interior point o (needed for a tetrahedronization 
of a convex polytope), after deleting points, lies inside the new convex hull, then the same 
formulas and time complexity, as by adding points, follow. If o lie outside the new convex 
hull, then, we need to choose a new interior point and recompute the new tetrahedra 
associated with it. Thus, we need in total 0{n) time to update the principal components. 

Under certain assumptions, we can recompute the new principal components faster: 
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• If we know that a certain point of the polyhedron will never be deleted, we can choose 
o to be that point. In that case, we also have the same closed-formed solution as 
for adding a point. 

• Let the facets of the convex polyhedron have similar (uniformly distributed) area. 
We choose o to be the center of gravity of the polyhedron. Then, we can expect 
that after deleting a point, o will remain in the new convex hull. However, after 
several deletion, o could lie outside the convex hull, and then we need to recompute 
it and the associate tetrahedra with it. 

Note, that in the case when we consider boundary of a convex polyhedron (Subsec- 
tion [6]T21 and Subsection I6.2.2p . we do not need an interior point o and the same time 
complexity holds for both adding and deleting points. 

6.1.2 Continuous PCA over a boundary of a polyhedron 

Let X be a polyhedron in M"^. We assume that the boundary of X is triangulated (if it 
is not, we can triangulate it in preprocessing), containing n triangles. The /c-th triangle, 
with vertices xi^k,X2,k,X3^k, can be represented in a parametric form by Tk{s,t) = xi^k + 
s {x2,k — xi^k) + 1 {x-s^k — xi^k), for < s, t < 1, and s + 1 < 1. For 1 < z < 3, we denote 
by Xij^k the i-th coordinate of the vertex Xj of the triangle T^. 
The center of gravity of the fc-th triangle is 



The contribution of each triangle to the center of gravity of the triangulated surface is 
proportional to its The area of the k-th triangle is 



We introduce a weight to each triangle that is proportional with its area, define as 




ttk = area(Tfc) 



\{x2,k - Xi^k)\ X |(x3,fc - fi,fc)| 



2 



ttk ttk 



En ! 
,=1 ak a 



where a is the area of X. Then, the center of gravity of the boundary of X is 



n 




k=l 



The covariance matrix of the k-th triangle is 



/o ft" iTk{s,t)-il) {n{s,t)-^l)T dtds 
Jo So'" 
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The {i,j)-th element of S^, i,j G {1,2,3}, is 

with fl = (/ii, /i2, /^s)- Finally, the covariance matrix of the boundary of X is 
Adding points 

We add points to X. Let X' be the new convex hull. We assume that X' is obtained 
from X by deleting rid, and adding Ua tetrahedra. Then the sum of the areas of all 
triangles is 

n na rid ria 

fc=l k=l k=l k=l k=l 

The center of gravity of X' is 

^^ = Efc=i w^fc/^fc + Efc=i w^fc/^fc - Ea,.=i ^fe/^fc 

= ^7 (ELi + Efcli Ofc/ifc - Efc=i Ofc/ife) (30) 

= ^7 (a/i + Efcli (^kflk - Efc=i akflk) ■ 

Let 

ria rid 

/ia = — Ofc/ifc, and /id = — flfc/ifc- 



a' ^ a' 

fc=l k=l 



Then, we can rewrite fl5U]) as 

fl' = —11 + fla- fid- (31) 
a' 



The z-th component of and fid, 1 < z < 3, is denoted by /ij_a and /ij^d, respectively. 
The (z, j)-th component, a-^, 1 < j < 3, of the covariance matrix S' of X' is 
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ii ( ELi E?=i ELi K{xi,i,k - f^'i){xj^h,k - /ij)+ 
ELi E?=i Ki^i,i,k - f^'i){xj,i,k - ^j] 

T2 ( Efcli E?=i ELi Ki^iAk - f^di^jAk - /ij)+ 
Efcli ELi Ki^i,i,k - fj''i)ixj,i,k - fj^'j)- 
Efc=i E?=i J2LiKixi,i,k - i^'i){xjxk - /ij)- 
Efeli E?=i Ki.Xi,i,k - f^'i)i^j,i,k - ^j] 

22 



Let 



where, 



/ ^ ( ' i' i' i' ' ' \ 



n 3 3 



= 5Z 5Z 5Z ^'ki.Xi,i,k - f^'i){Hh,k - /ip, (32) 

k=l 1=1 h=l 
n 3 

^i3,i2 = ^Y1 ~ f^'i)iXj,l,k -fJ''j), (33) 

k=l 1=1 

Ha 3 3 

^ij,21 = 5Z 5Z 5Z ^fc(^i.'.^' ~ ^^'^)^^hh,k - fi'j), (34) 
fc=l i=l /i=l 

ria 3 

(^ij,22 = ^Y1 ^fc(^i.'.fc ~ fJ''i)i^j,l,k -fJ''j), (35) 

fc=l 1=1 
"d 3 3 

^ii,3i = 5Z 5Z 5Z ^fc(^i.'.fe ~ f^'i)i^jAk - f^'j), (36) 
k=l 1=1 h=l 

rid 3 

^ij,32 = ^Y1 ^ki^i'^^k - fJ''i){Xj,i,k -fi'j)- (37) 
k=l 1=1 

Plugging-in the values of /i^ and fi'j in (132|) . we obtain: 

(^'ij,ll = ELl ELi ELi KiXi,l,k - ^f^i - f^i,a + f^i,d){Xj,h,k - ^f^j - f^j,a + f^j,d) 
= ELl E?=l ELi WkiXi,l,k - l^k + /^i(l - ^) - + f^i,d) 

{Xj,h,k - yUj + yUi(l - ^) - /^i,a + yUi,d) 

= ELl E?=i ELi <(a;i,;,fc - f^k)ixj^h,k -f^j)+ 

ELl E?=i ELi Kixi,i,k - /^i)(/^i(i - ^) - /^j,a + f^j,d)+ 
ELl E?=i ELi ^fc(/^i(i - ^) - + f^i,d)i^jAk - h)+ 

ELl E?=l ELi ^fc(/^i(l - ^) - f^i,a + /ii,rf)(/^i(l - ^) - /^i,a + 

(38) 

Since XlLi Ef=i ^'ki^i,i,k - fJ'i) = 0, 1 < i < 3, we have 
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- ^ J2l=l Ya=1 Yfh=l (^kiXi^i^k - fJ'i){Xj,h,k - fJ'j) + 

^ ELl ELi ELi - ^) - l^i,a + - ^) - f^j,a + 

= ELl Ef=i ELi (^k{xi,i,k - fj'i){xj,h,k - 

(39) 

Plugging-in the values of /i^ and /x^- in (155]) . we obtain: 

'^ij,12 = ELl ELi - f^f^i - f^i,a + f^i,d)i^j,h,k " " /^i,a + H,d) 

= ELl E?=l ^fc(a;i,«,fc - /^i + /^i(l - ^) - /^i,a + /ii,d) 
(a;j,h,fc - /ij + /ij(l - ^) - /ij,a + /ii,d) 

= ELiE?=i«^fcKi,fc (40) 

ELl E?=l W^fc(a^i,«,fc - - ^) - /^i,a + /ii,d) + 

ELl E?=l W'kifliil - ^) - /^i,a + fiLd)iXjM.,k - + 
ELl E?=l - ^) - + /^i,d)(/^i(l - ^) - + fij,d)- 

Since ELl ELi ^ki^i,i,k -/ii) = 0, l<z<3, we have 

'^ii,12 = ELl Ef=l 0'k{Xi,l,k - f^i) {Xj,h,k - fJ'j) + 

ELl Ef=l - ^) - /^i,a + fJ'i,d)it^ji^ ~ ^) ~ /^J> + , , 

3 *^ 

= ^ ELl E/=i o.kixi^i,k - /^i) (a^i,/i,fc - ^i)+ 

From (|40D and (|4T1). we obtain 



^ii,l = <ll+^ii,12 



(42) 

= aij + 12^(/ii(l - ^) - /ii,a + fli,d){flj{l - ^) - yUj-a + flj,d)- 

Note that cr^'^- j^ can be computed in 0(1) time. The components o"j'j- 2i and cr^ 22 can 
be computed in 0{na) time, while 0{nd) time is needed for computing Cj'j 31 and cr^j^^2- 
Thus, fj,' and 

~ ]^('^ij,ll + '^ij,12 + ^ij,21 + '^ij,22 + ^ij,31 + '^ii,32) 

= n(0-ij + a'ij 21 + 0-ij,22 + f^ij,31 + (^ij,32)+ (43) 
- ^) - f^i,a + /^i,d)(/^i(l - ^) - /^i,a + f^j,d)- 
24 



can be computed in 0{na + rid) time. 
Deleting points 

Let the new convex hull be obtained by deleting Ud tetrahedra from and added Ua 
tetrahedra to the old convex hull. Consequently, the same formulas and time complexity, 
as by adding points, follow. 



6.2 Continuous PCA in 

6.2.1 Continuous PCA over a polygon 

We assume that the polygon X is triangulated (if it is not, we can triangulate it in prepro- 
cessing), and the number of triangles is n. The k-th triangle, with vertices Xi^^, X2,fc, Xs,^ = 
o, can be represented in a parametric form by Tj(s, t) = x^ ^ + s (xi^^ — X3^fc)+t (x2,fc— a^s.fc), 
for < s,t < 1, and s + t < 1. 

The center of gravity of the fc-th triangle is 

'1 fl-s , 



lo lo ' ^^^^ 



The contribution of each triangle to the center of gravity of X is proportional to its area. 
The of the 2-th triangle is 

^rj. X I {X2,k - Xi^k) I X I (f3,fc - f I 

Ofc = area(Tfc) = 



2 

where x denotes the vector product. We introduce a weight to each triangle that is 
proportional with its area, define as 

Wk = ^ = — , 

where a is the area of X.Then, the center of gravity of X is 

n 

= ^ Wkflk- 
k=l 

The covariance matrix of the k-th triangle is 

y ^ /o ft" iTk{s,t)-^l) {n{s,t)~^l)T dtds 

Ej=i(^j,fc - /^)(^i,fc - l^f)- 

The (z, j)-th element of S^, i,j G {1, 2}, is 
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with fl = (/ii, /i2). The covariance matrix of X is 



Adding points 

We add points to X. Let X' be the new convex hulL We assume that X' is obtained 
from X by deleting n^, and adding Ua triangles. Then the sum of the areas of all triangles 
is 

k=l k=l k=l k=l k=l 

The center of gravity of X' is 

= Efc=i y^fc^fc + Efc=i w^fc/^fc - E A ^fe/^fc 

= ^7 (ELi + Efcli - Efc=i Ofc/ifc) (44) 

= (a/i + Efcli «fc/^A,. - Efc=i «fc/^fc) • 



Let 



f^a = — y~] akf^k, and y^d = — a-kf^k- 

k=l k=l 

Then, we can rewrite (144|) as 

/i' = -^yu + /ia - /id- (45) 
a' 

The 2-th component of and fid, 1 < i < 2, is denoted by /Xj ^ and /ij^d, respectively. 
The (z, j)-th component, o"-^, 1 < z, j < 2, of the covariance matrix E' of X' is 
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ii ( ELi ELi ELi Ki^iAk - f^'i)ixj,h,k - /ij)+ 

ELi E?=i Ki^i,i,k - f^'i)i^j,i,k - + 

T2 ( Efcli E?=i ELi Ki^iAk - f^di^jAk - /i^)+ 
Efcli ELi Wk{xi^i,k - fj''i)ixj,i,k - f^'j)- 
Er=i E?=i Y^i=iKixi,i,k - i^i){xjxk - /ij)- 

Efc=l E?=l Ki.^^,l,k - l^'di.^j,l,k - fJ'j, 



Let 
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where, 



n. 3 3 



^ij,n = XI XI 5Z '^ki^i'i''^ - f^'i)i^jAk - f^'j), (46) 

k=l 1=1 h=l 
n 3 

^ii,12 = X X ^'ki^i'l'k - fJ'diXj^i^k -fi'j), (47) 
k=l 1=1 

"a 3 3 

(^ij,2i = X X X ^fc(^^'.^ ~ /^')(3;i,/^,fc - f^'j), (48) 

fc=l i=l h=l 
ria 3 

f^ii,22 = X X ^fc(^^.'.fc ~ l^'i)i.^i,i,k - /^j), (49) 
fc=l i=l 

"d 3 3 

^ii,3i = X X X ^A,.(^^'.fc ~ f^'i)i^jAk - /ij), (50) 



fc=l i=l /i=l 
»id 3 



^ij,32 = X X ^fc(^^.^.fc ~ ^^'i}i^j,l,k - (51) 

fc=l 1=1 

Plugging- in the values of /i^ and n'^ in (H^ . we obtain: 

^ij,ll ~ X]fc=l X]«=l X]/i=l '^ki-'^i,l,k ~ '^f^i ~ f^i,a + f^i,d){Xj,h,k ~ '^f^j ~ f^j,a + A'-i.d) 
= ELl E?=l ELi ^fc(a;i,i,fc - + /^i(l - ^) - /^*,a + /i*,d) 

(a;j>,fc - /ij + /ij(l - f ) - fij,a + /ii,d) 

= ELl E?=l ELi KiXi,l,k - f^i)iXj,h,k - f^j) + 

ELl E?=l ELi 'WkiXi,l,k - /^i)(/^j(l - ^) - /^j,a + f^j,d) + 

ELl ELi ELi Kif^ii^ - ^) - + /^i,rf)(^i,/^,fc - n)+ 

ELl E?=l ELi «^fc(/^i(l - ^) - /^i,a + /^i,d)(/^j(l - ^) - /^i,a + /ii,d)- 

(52) 

Since Y2=i ELi w'k{^i,i,k -/ii) = 0, 1<2<2, we have 

"^ij.ii = ELl Ef=i ELi (^k{xi,i,k - fj'i){xj,h,k - 

^ ELl E?=l ELi «fc(/^i(l - ^) - f^i,a + /ii,d)(/ij(l - ^) - f^j,a + f^j,d) 

= ELl Ef=i ELi (^k{xi,i,k - fj'i){xj,h,k - fj'j)+ 

(53) 

Plugging- in the values of //^ and /^^ in fH7|) . we obtain: 
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= ELl E?=l KiXi,l,k - /ii + - f ) - /ii,a + 

= ELiE?=i«^fcK«,fc (54) 

ELl Ef=l W^fc(/^i(l - ^) - f^i,a + f^Ld)i^j,h,k - + 
ELl E?=l - ^) - /^i,a + - ^) - /^i,a + 

Since Y2=i ELi ^f^U^^i,;,^ -/ii) = 0,l<2<2, we have 



^ij,12 = ^7 ELl Ef=l (^k{Xi,l,k - P-i) {Xj,h,k - fJ'j) + 

^ ELl ELi - ^) - ^i,a + - ^) - /^i,a + /^i.d) 

= ^ ELl Ef=i 0'k{xi,i,k - f^i){xj,h,k - fj'j)+ 

From (1541) and (ISSj) . we obtain 



(55) 



(56) 

= fly + 12^(/i,(l - f ) - /i.,a + /i*,d)(yUj(l - ^) - /^i,a + f^j,d)- 

Note that cr^'^ ^ can be computed in 0(1) time. The components cr^'^ gi and cr-j^22 can 
be computed in 0(na) time, while 0{nd) time is needed for computing cr'ij^i and cr^j^^2- 
Thus, fj,' and 



= ]I«,ll+<,12 + <,21+<,22 + <,31 + <, 



= + ^ij,21 + <-,22 + f^ij,31 + (^ij,32)+ (57) 

can be computed in 0{na + Ud) time. 
Deleting points 

Let the new convex hull be obtained by deleting Ud tetrahedra from and added Ua 
tetrahedra to the old convex hull. If the interior point o, after deleting points, lies inside 
the new convex hull, then the same formulas and time complexity, as by adding points, 
follow. However, o could lie outside the new convex hull. Then, we need to choose a new 
interior point o', and recompute the new tetrahedra associated with it. Thus, we need in 
total 0{n) time to update the principal components. 
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6.2.2 Continuous PC A over the boundary of a polygon 

Let X be a polygon in M^. We assume that the boundary of X is comprised of n line seg- 
ments. The k-th line segment, with vertices Xi^k, 2:2, fc, can be represented in a parametric 
form by 

Skit) = Xi^k + t ix2,k - 

Since we assume that the mass density is constant, the center of gravity of the k-th line 
segment is 

/q Skit) dt _ xi^k + X2,k 



lo dt 



10 

The contribution of each line segment to the center of gravity of the boundary of a polygon 
is proportional with the length of the line segment. The length of the k-th line segment 
is 

Sk = length(S'fc) = p2,fc - 
We introduce a weight to each line segment that is proportional with its length, define as 

Sk Sk 

Wk = = — , 

where s is the perimeter of X. Then, the center of gravity of the boundary of X is 

n 

= y^^wkfik- 

k=l 

The covariance matrix of the fc-th line segment is 

^k — (-1 jj 

= I ( E?=i ELi i^hk - P) i^h,k - 

The (i,j)-th element of S^, i.j G {1,2}, is 

^ij,k = g ^ E«=l E/i=l("^*>'>fc ~ f'i)i-^j,h,k ~ f^j)~ 
YM=li^i,l,k — ^■i)iXj,l,k — fJ'j) 

with fl = (/ii,/i2)- 

The covariance matrix of the boundary of X is 
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Adding points 

We add points to X. Let X' be the new convex hull. We assume that X' is obtained 
from X by deleting Ud, and adding Ua line segments. Then the sum of the lengths of all 
line segments is 

n ria Ua 

fc=l fc=l k=l k=l k=l 

The center of gravity of X' is 

= 7 (ELi ^fc/^fc + Efcli 5fc/^fc - Efc=i ^fc/^fc) (58) 

Let 

yWa = — ^ SfcyUfc, and f^d = — ^ Sheik- 



s' ^ s' 

k=l k=l 



Then, we can rewrite fISS]) as 



fl' = + fla- fid- (59) 



The i-th component of fla and fid, 1 < z < 2, is denoted by /ij^a and /ij^^, respectively. 
The (z,j)-th component, a-^, 1 < j < 2, of the covariance matrix S' of X' is 

(^'ij = l{ ELi ELi ELi Kixi,i,k - /iD (a;j,h,fc - + 
ELi ELi w^fc(a;i,«,fc - f^dixj,i,k - /x;)) + 

i(Efcli ELi J2l=iKixi,i,k - fj''i){xj,h,k - 
TIU ELi ELi K{xi,i,k - f^dixj,h,k - /ij)- 

Efc=i E?=i Wfc(a;i,«,fc - f^'i)ixj,i,k - /ij)) . 

Let 



where, 



2 2 



= XI XI XI ^fc(^^'.fc ~ f^'i)i^jAk - fJ^'j), (60) 



fc=i 1=1 h=i 
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n 2 

k=l 1=1 

Ua 2 2 

^'i3,2\ = Y^Y^^ W'ki^i,l,k - f^'i)iXj,h,k - /ij), 
k=l 1=1 h=l 

ria 2 

(^'ij,22 = XI XI ^'ki^i'l'k - ^^'i}{Xj,l,k - /^j), 
k=l 1=1 

na 2 2 

= X X X ~ f^'i)(^jAk - fJ^'j), 

k=l 1=1 h=l 

rid 3 

'^ij,32 = ^^1^kiXi,l,k - f^'i)i^j,l,k - tJ'j)- 
k=l 1=1 

Plugging-in the values of /i^ and fi'j in fl60|) . we obtain: 

l^i){Xj,h,k - + 

Since XlLi SLi ^'ki^i,i,k - Hi) = 0, 1 < i < 2, we have 
^ii.ii = 7 ELi E/=i ELi Sfc(a;i,i,fc - fii){xj,h,k - + 

7 ELl E/=l ELi Sfc(/ii(l - 7) - f^i,a + - 7) - /^i,a + fij,d) 

= 7 ELi E?=i ELi - ^i^){xj,h,k - 

Plugging-in the values of /i^ and /x^- in (16T]) . we obtain: 





.EL 


=1 E/t= 




2^k= 


.EL 


=1 E/i= 










iS^j,h,k 


l^k= 


.EL 


=1 Eh= 


-_l Wk{Xi^l,k - 


s--\n 
l^k= 


.EL 


=1 Eh= 


-.l'^ki^i,l,k - 


l^k= 


.EL 


=1 Eh,= 


-.1 - 


\--\n 
l^k= 


.EL 


=1 E/i= 
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{Xj,h,k - fJ'j + - 7) - fJ'j,a + 

= 'Zt=iJ2LiKi^i,i,k- f^i){xj,h,k- (68) 

ELl E?=i^fe(/^i(i - 7) - /^i,'^ + /^i,<i)(/^i(i - 7) - /^J> + 

Since ELi ELi '^'ki.^hhk -/ij) = 0, 1<?<2, we have 



(^'ij,l2 = 7 ELl ELi Sk{Xi^l,k - /ii) - + 

7ELiELi*fc(/^i(l-7)-/^i." + /^M)(/^i(l-7)-/^J> + /^i.<^) , , 
, (69) 

= 7 ELl ELi Sk{xi,i,k - f^i){xj^h,k - ^i)+ 

2Jr(/ii(l - 7) - /ii,a + /ii,d)(/ij(l - 7) - /ij,a + /^i.d)- 

From (EHD and (l69l). we obtain 



(70) 

= CTjj + 6pr(/-ii(l - 7) - /ii,a + Ati,d)(/ii(l - 7) - Atj,a + /^i.d)- 

Note that jl' and the components o'[j2i^ ^ij, 22^ ^i,3i '^'1,32 '^^^ ^e computed in 
constant time, under the assumption that X and X' differ in the constant number of 
polyhedra, i.e., Ua and Ud are constants. In that case, also the component 

'^ij ~ ^i'^ij^l ~^ '^ijA2 ~^ '^ij,21 ~^ '^ij,22 ~^ '^ij,31 ~^ '^ij,32) 

= lia^j + a^21 + C^ij,22 + <-,31 + ^ii,32)+ (71) 

jrifiiil - jr) - f^i,a + /ii,d)(/ii(l - 7) - f^j,a + /^i,d) 
can be computed in 0{na + Ud) time. 

Deleting points 

Let the new convex hull be obtained by deleting Ud tetrahedra from and added 
tetrahedra to the old convex hull. Consequently, the same formulas and time complexity, 
as by adding points, follow. 
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